#define NTHREAD 20
#define NGR		105

__constant__ unsigned int shift1[4] = {6, 2, 13, 3};
__constant__ unsigned int shift2[4] = {13, 27, 21, 12};
__constant__ unsigned int shift3[4] = {18, 2, 7, 13};
__constant__ unsigned int offset[4] = {4294967294, 4294967288, 4294967280, 4294967168};

__shared__ unsigned int randStates[32];

__device__ unsigned int TausStep(unsigned int &z, int S1, int S2, int S3, unsigned int M)
{
    unsigned int b = (((z << S1) ^ z) >> S2);
    return z = (((z &M) << S3) ^ b);
}

__device__ unsigned int randInt()
{
    TausStep(randStates[threadIdx.x&31], shift1[threadIdx.x&3], shift2[threadIdx.x&3],shift3[threadIdx.x&3],offset[threadIdx.x&3]);
    return (randStates[(threadIdx.x)&31]^randStates[(threadIdx.x+1)&31]^randStates[(threadIdx.x+2)&31]^randStates[(threadIdx.x+3)&31]);
}


extern "C" __global__ void perturbSolution(float** graph, int* currentSol, int** solutionSpace, int numPerturb, 
										   int n, float *t, int *stage, int* flags)
{
   
   int tx = threadIdx.x;
   int tidX = blockDim.x * blockIdx.x + threadIdx.x;
   if(tidX>=256)
     return;
   int fPtr = tidX*n;
   __shared__ int path[NTHREAD][NGR];
   __shared__ int mystage;
   mystage = stage[0];
   int pack = -1;
   float COST2;
   float COST1=0;
   for(int i=0;i<n;++i) {

      path[tx][i] = currentSol[i];
   }   
   
   
   for(int i=0;i<n;++i) {
   	  COST1 += graph[path[tx][i]][path[tx][(i+1)%n]];
   	  
   }
   
   for(int ix = 0;ix<numPerturb; ++ix) {
   	  if(stage[0] == 1) {
   	  	for(int i=0;i<tidX/32;++i) {
   	  		if(flags[i]!=0)
   	  			goto end;
   	  	}
   	  }
      unsigned int k = randInt()%n;
      unsigned int kk = randInt()%n;
      int TEMP = path[tx][k];
      path[tx][k] = path[tx][kk];
      path[tx][kk] = TEMP;
      
      COST2 = 0;
      
      for(int j=0;j<n;++j) {
   	  	COST2 += graph[path[tx][j]][path[tx][(j+1)%n]];
   	  	
      }
	 
      if(COST2<COST1) {
      		
      		if(stage[0]==1) {
	      		pack = tidX/32;
	      		int offset = 32-tidX%32;
	      		int mask = 1;
	      		mask = mask <<(offset-1);
	      		flags[pack] = flags[pack]|mask;
	      	}
	      	COST1=COST2;
      }
      else
      {
      	 
      	  float prob =  ((randInt()%10000))/10000.0;
      	  float delta = COST2-COST1;
      	  if(prob<(exp(-delta/t[0]))) {
      		COST1=COST2; 
      		if(stage[0]==1) {
	      		pack = tidX/32;
	      		int offset = 32-tidX%32;
	      		int mask = 1;
	      		mask = mask<<(offset-1);
	      		flags[pack] = flags[pack] | mask;
	      	}
      	  }
      	  else
      	  {
      	  	TEMP  = path[tx][k];
      		path[tx][k] = path[tx][kk];
      		path[tx][kk] = TEMP;
      	  }
      } 
      if(stage[0]==1 && pack!=-1) {
      	goto end;
	  }
   }
   end:
   for(int z=0;z<n;++z) {
    	solutionSpace[tidX][z] = path[tx][z];
   }
}

	